# Name:         IceModeling.py script 
#
# Content:      Relevant functions for ice modeling
#
# Author:       Ben Martinez
# E-mail:       benjamin.martinez@vattenfall.com
# Company:      Vattenfall R&D, Fredericia, DK
# Copyright:    (c) bejm 2012
# Licence:      <vattenfall>
#
#-------------------------------------------------------------------------------

# Adding path to windroseicing script & loading modules

import sys

Path_windroseicing = 'C:\Ben\R&D Documentation\icing\code\windrose'
sys.path.append(Path_windroseicing) 

import numpy as np
import math
from scipy.ndimage import map_coordinates
import matplotlib.pyplot as plt
from windroseicing import WindroseAxes
from matplotlib import mpl
import datetime as dti

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

def IceModel(Vmag,w,T,p,rh,Nc,dt,Dstart):

# Function:     IceModel
#
# Purpose:      Ice model providing cylinder-based ice load estimation [kg/m].
#               To be used with standard mesoscale model outputs as input.
#
# Assumptions:  Freely rotating cylinder, ice accretion evenly distributed,
#               considers dry growth and density variation of the accreted ice.
#           
# Created:      21/12/2012
#
#------------------------------ USAGE ------------------------------------------
#
#    IL = IceModel(Vmag,w,T,p,rh,Nc,dt,Dstart)
#
#--- input arguments (in the form of numpy arrays):
#
# Vmag = Horizontal wind speed magnitude (x,y components) time series array [m/s]
# w = Cloud liquid water content time series array in [g/kg] (per kg of dry air)
# T = Temperature time series array [Celcius]
# p = Atmospheric pressure time series array [hPa]
# rh = Relative humidity time series array in [%]
# Nc = Cloud droplet concentration time series array in [cm^-3]
# dt = Length of time step [s] , constant
# Dstart = Initial diameter of cylinder [m] (normally 30*1e-3 [m]), constant  
#
#--- output (in the form of a list) 
#
# IL = ice load time series in [kg/m]
#
#----------------------------- REFERENCES --------------------------------------
#
# [1] Finstad et alt. 1988, J Oceanic Atmos. Technol., 5, 160-170
# [2] ISO 12494, Atmospheric icing on structures
# [3] Roland B. Stull, Meteorology for Scientists and Engineers
# [4] Bernstein B. et alt, Personal comunication
# [5] Makkonen L., Modeling of Ice Accretion on Wires
# [6] Makkonen L., Models for growth of rime, glaze, icicle and 
#     wet snow on structures
# [7] Harstveit K., Using Routine Meteorological data from airfields 
#     to produce a map risk zones in Norway
#
#---------------------------- Constants -------------------------------------

    Rd = 287.04  # Specific dry air gas constant [J.K^-1.kg^-1]    
    Kel = 273.15  # Degrees Celcius to Kelvin offset [K] 
    rho_w = 1000.  # Water density [kg/m^3]
    alpha1_min1 = 0.01 # collision efficiency lower limit [1]
    alpha1_min2= 0.005 # Absolute lower limit 
    alpha_2 = 1. # Sticking efficiency of water droplets approximated to 1.
                 # (Non-bouncing droplets is a good aprox. [2])             
    e_o = 611. # Tetens formula: reference water-vapor saturation pressure [Pa]
    bs = 17.2694 # Tentens formula: b parameter [-]
    T2s = 35.86 # Tetens formula: reference temperature [K]
    ka = 0.023   # Thermal conductivity of air [W/(m.K)]
    Lf = 334000.  # Latent heat of (fusion/melting) [J/kg]
    lamb = 0.26   # Liquid fraction during wet growth [-]
    Ts = 273.15  # Temperature of the icing surface in wet growth [K]
    eps = 0.622  # ratio of molecular masses of dry air and water vapour [-]
    Le = 2500000.  # Latent heat of evaporation/condensation [J/kg]
    cp = 1004.67  # Specific heat of dry air [J/(kg.K)]
    cw = 4200.    # Specific heat capacity of water [J/(kg.K)]
    e_s_surf = 611 # Saturation water vapour pressure over accretion surface(0 Cel)[Pa]
    r = 0.79    # Recovery factor for viscous heating on a cylinder
    sigma = 5.67*1e-8    # Stefan-Boltzmann constant [W/(m^2.K^4)]
    a_r = 8.1*1e7   # Radiation linearization constant [K^3]
    subl_factor = 0.02  # Sublimation factor (originally 0.2) [-]
    shed = 4.   # Emperical Shedding factor [-] (To be tunned)

#---------------------------- Initialization --------------------------------#

    jt = len(T)  # Time series length
    Mice = 0     # Initial total ice mass [kg/m]
    rho_it = 0   # Initial estimated density of total ice mass [kg/m^3]
    Dm = Dstart  # Inicial droplet diameter [m]
    ice_load = [Mice]  # Inicial ice load list [kg/m]

#---------------------------------------- -----------------------------------#

    for j, i in enumerate(T):

        # Pressure in Pascals [Pa]       
        pr = p[j]*100.  

        # Air density calculation:
        rho_a = pr/((T[j] + Kel)*Rd)  # Atmospheric air density [kg/m^3]

        # Converting LWC to kg/m^3 units 
        #LWC  = w[j]*rho_a*1e-3  #(input given in g/kg of dry air)
        LWC = w[j]*1e-3 # (select this if input given in g/m^3)

        # dynamic viscosity of dry air temperature dependance
        mu_a = 1.718*1e-5 + 5.1*1e-8*T[j]
     
#------------------- Computing Collision efficiency alpha_1 ------------------#    
    
        if w[j] > 0. :
        
            mu_s = min(np.array([1000./Nc[j] + 2., 15]))
    
            lambdas = ((np.pi/6.)*rho_w*(math.gamma(4. + mu_s)/math.gamma(1. + mu_s))*\
                  (Nc[j]*1e6/LWC))**(1./3.)

            # Median Volume diameter in [m]
    
            MVD = (3.672 + mu_s)/lambdas
          
            #Droplet Reynolds number
    
            Re = (rho_a*MVD*Vmag[j])/mu_a
        
            # K and Phi dimensionless parameters  
        
            K = (rho_w*Vmag[j]*(MVD**2.))/(9*mu_a*Dm)
            phi = (Re**2.)/K
        
            # A, B and C colision efficiency calculation
        
            A = 1.066*(K**(-0.00616))*np.exp(-1.103*(K**(-0.688)))
            B = 3.641*(K**(-0.498))*np.exp(-1.497*(K**(-0.694)))
        
            if phi > 100.:
                C = 0.00637*((phi-100.)**(0.381))
            else:
                C = 0.
            
            # Collision efficiency alpha_1:
        
            alpha_1 = A - 0.028 - (C*(B - 0.0454))
        
            if alpha_1 < alpha1_min1:
                alpha_1 = alpha1_min2/(Dm**(0.5))
            
        else:
         
            alpha_1 = 0.
            
#-------------------- Computing Accretion efficiency alpha_3 ----------------#
            
        # Saturation pressure of water-vapor (Tetens formula, [3])
    
        e_s = e_o*np.exp((bs*T[j])/(T[j] + Kel - T2s))
            
#--- Estimating the convective heat transfer coefficient (h) of a cylinder---#

        # Reynolds number of the cylinder
    
        Re_c = (rho_a*Dm*Vmag[j])/mu_a    

        # Nusselt number [7]   
    
        Nu = 0.032*(Re_c**(0.85))

        # Convective heat transfer coef. [7]
    
        h = (ka*Nu)/Dm
    
    
        if w[j] > 0. :
        
            if T[j] <= 0. and T[j] > -4. :
        
            # Flux density of water to surface
    
                F = alpha_1*alpha_2*LWC*Vmag[j]
        
                alpha_3 = (h/(F*Lf*(1.-lamb)))*( (Ts - (T[j] + Kel)) + \
                       (((eps*Le)/(cp*pr))*(e_s_surf - e_s)) - \
                       ((r*(Vmag[j]**2))/(2.*cp)) ) + \
                       ((cw*(Ts - (T[j] + Kel)))/((1.-lamb)*Lf)) + \
                       ((sigma*a_r*(Ts - (T[j] + Kel)))/(F*Lf*(1.-lamb)))
    
                if alpha_3 > 1. :        
                    alpha_3 = 1.
            
                elif alpha_3 < 0 :            
                    alpha_3 = 0.
          
            elif T[j] <= -4. :
                alpha_3 = 1.
            
            else:
                alpha_3 = 0.
            
          
        else:
            alpha_3 = 0.
            
#---------------- Computing the ice mass for each time step ------------------#

        # - Computing accreted mass factor only - #      
    
        if T[j] < 0. :     
            dMdt = alpha_1*alpha_2*alpha_3*LWC*Vmag[j]*Dm
    
        else:
            dMdt = 0.
          
        # - Computing Sublimation based on Ben Bernstein et alt [4] - #
    
        if rh[j] >= 90. :
            rhmap = 1.
        elif rh[j] >= 80. and  rh[j] < 90. :
            rhmap = 0.25 + (0.75*(rh[j] - 80) / 10.)
        elif rh[j] >= 70. and rh[j] < 80. :
            rhmap = 0.25*(rh[j] - 70.)/10.         
        else:
            rhmap = 0.
        
        if dMdt > 0. :
            subl = 0.
        
        else:
            subl = subl_factor*(dt/3600.)*((0.65*(min(1.,Vmag[j]/10.))) + \
                  (0.35*(1.-rhmap)))
              
        # - Computing Melting and Shedding factors all together - #          
              
        if T[j] > 0. and  Mice > 1e-20 :       
            # Available Convective heat              
            Qh = h*T[j]*shed  
     
            # Available Evaporative heat       
            Qe = ((e_s - e_s_surf)*h*eps*Le*shed)/(cp*pr)
            
            if Qe < 0. :
                Qe = 0.
     
            # Available Radiation heat      
            Qr = sigma*shed*(((Kel + T[j])**4.) - Kel**4.)
     
            # Total heat for melting     
            Q_melt = Qh + Qe + Qr
        
            # Melted mass of ice (assuming Q_melt = Q_fusion)
            melt = (Q_melt*np.pi*Dm*dt)/Lf
        
        else:
            melt = 0.
            
        # - Computing total accreted ice and updating diameters and densities - # 

        if j < (jt-1) : 

            Delta_Mice =  Mice + dMdt*dt - subl - melt       
            Mice_new = max(0., Delta_Mice)
                
            # - Melting/Sublimation condition - #
            if (Mice_new < Mice) and (Mice_new > 1e-20) :
            
                # Update cylinder diameter and density
                Dm_new = (((4.*(Mice_new - Mice))/(np.pi*rho_it)) + \
                            Dm**2.)**(0.5)
                rho_it_new = rho_it            
            
            # - Accretion condition - #
            elif Mice_new > (Mice + 1e-15) :
            
                # Calculate new accreted ice density
                Ko = K*((0.087*(Re**(0.76*(Re**(-0.027)))) + 1.)**(-1.))
         
                if Ko > 0.55 :
                    Vo = Vmag[j]*(0.561 + 0.592*np.log10(Ko) - 0.26*((np.log10(Ko))**2.))
                else:
                    Vo = Vmag[j]*(-0.174 + 1.464*Ko - 0.816*(Ko**2.))
             
                T_s = T[j]  # Reasonable assumption according to Makonnen [5] 
         
                # Macklin's parameter
                R_ma = -1.*((Vo*(MVD*1e6))/(2.*T_s))
         
                if R_ma <= 0. :
                    R_ma = 0.0001
                
                # New accreted ice density in kg/m^3
                rho_ice_new = (0.378 + (0.425*(np.log10(R_ma))) - \
                              0.0823*((np.log10(R_ma))**2.))*1000. 
                              
                # Ice density lower limit according to Makkonnen [6]            
                if rho_ice_new < 100. :
                    rho_ice_new = 100.
                
                # Updating cylinder diameter using accreted ice density    
                Dm_new = (((4.*(Mice_new - Mice))/(np.pi*rho_ice_new)) + \
                            Dm**2.)**(0.5)  
                
                # Updating/averaging total ice density
                rho_it_new = (4.*Mice_new)/(np.pi*((Dm_new**2.)-(Dstart**2.)))                 
                            
            
            # - No mass on the cylinder - #
            elif Mice_new == 0. :
            
                Dm_new = Dstart
                rho_it_new = 0.
            
            # - Quasi-static conditions - #
            else:
        
                Dm_new = Dm
                rho_it_new = rho_it
            

            # - Creating ice load vector - #
            ice_load.append(Mice_new)

    
        # - Updating ice mass, cylinder diameter and ice density for new loop
        Mice = Mice_new
        Dm = Dm_new
        rho_it = rho_it_new
    
    return ice_load

#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

def IcePowerLossModel(Vmag,IL,PathPowerCurve,PathPowerCurve3D):
    
# Function:     IcePowerLossModel
#
# Purpose:      Calculates Aerodynamic Performance losses due to icing based
#               on a predefined 3D power curve. 
#
# Assumptions:  3D power curve is unique. Specific assumptions are made  
#               in order to accomodate all type of power curve shapes.
#               Within these, the most relevant:
#               - Percentage of performance loss is considered equal 
#                (input power curve with respect to 3D power curve)
#                for a given wind speed                              
#
# Created:      17/01/2013
#
#------------------------------ USAGE -----------------------------------------
#
#    function: 
#    Per_loss_list, P_with_ice_list, P_no_ice_list, P_loss_list = 
#    IcePowerLossModel(Vmag,IL,PathPowerCurve,PathPowerCurve3D)
#
#--- input arguments:
#
#   Vmag = Horizontal wind speed magnitude (x,y components) time series 
#         1D array [m/s]
#   IL = Ice load time series 1D array [kg/m]
#   PathPowerCurve = Path to Wind Turbine Power Curve file. Power Curve in 
#                  standard WAsP power curve format (.pow, ASCII) 
#   PathPowerCurve3D = Path to iced 3D power curve file. Stored in surfer 
#                      .grd format 
#
#--- output (in the form of lists) 
#
#   Per_loss_list = Percentage of power loss [%]
#   P_with_ice_list = Power output considering ice effects [kW]
#   P_no_ice_list = Power output without considering ice effects [kW]
#   P_loss_list = Total production loss due to icing [kW]

#------------------------- Reading input files --------------------------------

    # Reading standard WAsP power curves (.pow ASCII format)

    f = open(PathPowerCurve,'r')
    lines = f.readlines()
    
    line2 = lines[2].split()
    speed_scaling , Power_scaling = float(line2[0]) , float(line2[1])

    WS_pc_d = []
    PW_pc_d = []
    
    for j, i in enumerate(lines[3:]):
    
        line = i.split()
        WS_pc_d.append(float(line[0])*speed_scaling)
        PW_pc_d.append(float(line[1])*Power_scaling)
    
    WS_cut_off = WS_pc_d[-1] # Cut-off wind speed of the power curve
    WS_cut_in = WS_pc_d[0]   # Cut-in wind speed of the power curve
   
# ------------- Reading 3D iced power curve in .grd file ---------------------#

    g = open(PathPowerCurve3D,'rb')

    lines_g = [lin.split() for lin in g]
    xlines , ylines = int(lines_g[1][0]) , int(lines_g[1][1])
    WSmin , WSmax = float(lines_g[2][0]) , float(lines_g[2][1])
    ILmin , ILmax = float(lines_g[3][0]) , float(lines_g[3][1])
        
    Power_ice = np.array([map(float,lin) for lin in lines_g[5:]])

    # Other definitions

    WS_cut_in_ice = 3. # Cut-in wind speed for the iced power curve
    inc = 0.5 # Interval assumption when accomodating power curve with 
              # WS_cut_in < WS_cut_in_ice 

#---------------------------------------------------------------------------- #

    Per_loss_list = []
    P_with_ice_list = []
    P_no_ice_list = []
    P_loss_list = []
   
    for i , j in enumerate(Vmag):
    
        if (j > WS_cut_in) and  (WS_cut_in >= WS_cut_in_ice):

        
            if j >= WS_cut_off:
                P_no_ice = 0.              
            
            # Obtaining Power from power curve when IL = O kg/m
            else: 
                P_no_ice = np.interp(j,WS_pc_d,PW_pc_d)
            
                # Setting upper limit equal to WSmax (Assuming same loss when
                # WS > WSmax and WS_cut_off > WSmax
                if j > WSmax :
                    j = WSmax
                
            # Setting Ice load limits for interpolation
        
            if IL[i] < ILmin:
                IL[i] = ILmin
            
            if IL[i] > ILmax:
                IL[i] = ILmax
            
            # Converting to indices: x direction
            # (WS_int should range from 0 to (xlines - 1))
            WS_int = [(xlines - 1) * (j - WSmin) / (WSmax - WSmin), xlines - 1]

            # Converting to indices: y direction
            # (IL_int should range from 0 to (ylines - 1))
            IL_int = [(ylines - 1) * (IL[i] - ILmin) / (ILmax - ILmin), ylines - 1]

            if IL[i] > 0 :
            
                # Obtaining iced power ratio from iced power curve using map 
                # interpolation
                r_ice = map_coordinates(Power_ice, [IL_int, WS_int],order = 5)
                r_no_ice = map_coordinates(Power_ice,[[0,300] , WS_int], \
                order = 5)
                             
                # Percentage of power loss
                Per_loss  = (1 - (r_ice[0]/r_no_ice[0]))*100.
                
                # Power with ice
                P_with_ice = P_no_ice*(1 - (Per_loss/100.)) 
                
                # Power loss due to ice            
                P_loss = P_no_ice - P_with_ice

            else:
                P_loss = 0.
                Per_loss = 0.
                P_with_ice = P_no_ice
            
        elif (j > WS_cut_in) and (WS_cut_in < WS_cut_in_ice):
          
            if j >= WS_cut_in_ice:
            
                if j >= WS_cut_off:
                    P_no_ice = 0.              
            
                # Obtaining Power from power curve when IL = O kg/m
                else:
                    P_no_ice = np.interp(j,WS_pc_d,PW_pc_d)
            
                    # Setting upper limit equal to WSmax (Assuming same loss when
                    # WS > WSmax and WS_cut_off > WSmax
                    if j > WSmax :
                        j = WSmax
                
            # Setting Ice load limits for interpolation
        
                if IL[i] < ILmin:
                    IL[i] = ILmin
            
                if IL[i] > ILmax:
                    IL[i] = ILmax
            
                # Converting to indices: x direction
                # (WS_int should range from 0 to (xlines - 1))
                WS_int = [(xlines - 1) * (j - WSmin) / (WSmax - WSmin), \
                         xlines - 1]

                # Converting to indices: y direction
                # (IL_int should range from 0 to (ylines - 1))
                IL_int = [(ylines - 1) * (IL[i] - ILmin) / (ILmax - ILmin), \
                         ylines - 1]

                if IL[i] > 0 :
                
                    # Obtaining iced power ratio from iced power curve using 
                    # map interpolation
                    r_ice = map_coordinates(Power_ice, [IL_int, WS_int], \
                            order = 5)
                    r_no_ice = map_coordinates(Power_ice,[[0,300] , WS_int], \
                               order = 5)
                             
                    # Percentage of power loss
                    Per_loss  = (1 - (r_ice[0]/r_no_ice[0]))*100.
                    
                    # Power with ice
                    P_with_ice = P_no_ice*(1 - (Per_loss/100.)) 
                    
                    # Power loss due to ice            
                    P_loss = P_no_ice - P_with_ice

                else:
                    P_loss = 0.
                    Per_loss = 0.
                    P_with_ice = P_no_ice
            
            else:
            
                P_no_ice = np.interp(j,WS_pc_d,PW_pc_d)
            
                if IL[i] < ILmin:
                    IL[i] = ILmin
            
                if IL[i] > ILmax:
                    IL[i] = ILmax
                       
                WS_int = [(xlines - 1) * ((WS_cut_in_ice + inc) - WSmin) / \
                    (WSmax - WSmin), xlines - 1]
            
                IL_int = [(ylines - 1) * (IL[i] - ILmin) / (ILmax - ILmin), \
                         ylines - 1]
            
                if IL[i] > 0 :
                
                    # Obtaining iced power ratio from iced power curve using 
                    # map interpolation
                    r_ice = map_coordinates(Power_ice, [IL_int, WS_int], \
                            order = 5)
                    r_no_ice = map_coordinates(Power_ice,[[0,300] , WS_int], \
                               order = 5)
                             
                    # Percentage of power loss
                    Per_loss  = (1 - (r_ice[0]/r_no_ice[0]))*100.
                    
                    # Power with ice
                    P_with_ice = P_no_ice*(1 - (Per_loss/100.)) 
                    
                    # Power loss due to ice            
                    P_loss = P_no_ice - P_with_ice

                else:
                    P_loss = 0.
                    Per_loss = 0.
                    P_with_ice = P_no_ice
                        
        else:
        
            P_loss = 0.
            Per_loss = 0.
            P_no_ice = 0.
            P_with_ice = P_no_ice
       
        Per_loss_list.append(Per_loss)
        P_with_ice_list.append(P_with_ice/1000.)
        P_no_ice_list.append(P_no_ice/1000.)
        P_loss_list.append(P_loss/1000.)
        
    return [round(i,2) for i in Per_loss_list], \
    [round(i,2) for i in P_with_ice_list],[round(i,2) for i in P_no_ice_list] \
    ,[round(i,2) for i in P_loss_list] 
    
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

def IceRate(IL,dt):
    
# Function:     IceRate
# 
# Purpose:      Calculates Ice rate and Icing Active hours from input Ice 
#               load [kg/m] (as numpy array) and time step 
#               (constant given in s.). Result in g/(m.h)
#
# Created       17/01/2013

    
    IR = []
    for i , j in enumerate(IL[0:-1]):
        
        if IL[i + 1] > IL[i] :
            
            it = ((IL[i + 1]) - IL[i])/(dt/3600.)

        else:
            it = 0.
            
        IR.append(it*1000.)

    ActiveHours = len([i for i, e in enumerate(IR) if e >= 10.])
    
    return [np.round(k) for k in IR], ActiveHours
 
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# Plot functions

def  ActiveHoursVersusTemp(T,IR):
    
# Function: ActiveHoursVersusTemp   
#
# Purpose:  Plots Histogram of Temperature field and Active Icing Hours
#           Inputs are Temperature (T) and Ice Rate (IR)
#
# Created:  21/01/2013
    
    ind = []

    for i , j in enumerate(IR):
        if j >= 10:
            ind.append(i)
            
    T = np.array(T)
    wei1 = np.ones(len(T))*(100./len(T))
    wei2 = np.ones(len(T[ind]))*(100./len(T))
    hs = np.arange(int(min(T))-2,int(max(T))+1)
    
    plt.figure()
    plt.hist(T,hs,weights = wei1,color = 'g',label = 'Total time')
    plt.hist(T[ind],hs, weights = wei2,color = 'b',label = 'Active icing')
    plt.legend(loc = 'best')
    plt.xlabel('Temperature [deg]')
    plt.ylabel('Percentage of time [%]')
    plt.grid()     
    plt.show()
    
#-------------------------------------------------------------------------------    
    
def ActiveHoursVersusWindSpeed(WS,IR):
    
# Function: ActiveHoursVersusWindSpeed  
#
# Purpose:  Plots Histogram of Wind speed field and Active Icing Hours
#           Inputs are Wind Speed (WS) and Ice Rate (IR)
#
# Created:  21/01/2013    

    ind = []

    for i , j in enumerate(IR):
        if j >= 10:
            ind.append(i)    
            
    WS = np.array(WS)
    wei1 = np.ones(len(WS))*(100./len(WS))
    wei2 = np.ones(len(WS[ind]))*(100./len(WS))
    hs = np.arange(int(min(WS))-2,int(max(WS))+1)
    
    plt.figure() 
    plt.hist(WS,hs,weights = wei1,color = 'g',label = 'Total time')
    plt.hist(WS[ind],hs, weights = wei2,color = 'b',label = 'Active icing')
    plt.legend(loc = 'best')
    plt.xlabel('Wind Speed [m/s]')
    plt.ylabel('Percentage of time [%]')
    plt.grid()     
    plt.show()
            
#-------------------------------------------------------------------------------

def ActiveIcingRose(Wdir,IR):
    
# Function: ActiveIcingRose
#
# Purpose:  Plots Active icing rose in terms of probability of occurence
#           Inputs are Wind Direction (Wdir) and Ice Rate (IR) 
#
# Note:     To be used with len(IR) = len(Wdir) - 1 ; The code uses Wdir[1:]
#
# Created:  30/01/2013

    def new_axes():
        fig = plt.figure(figsize=(8, 8), dpi=80, facecolor='w', edgecolor='w')
        rect = [0.1, 0.1, 0.8, 0.8]
        ax = WindroseAxes(fig, rect, axisbg='w')
        fig.add_axes(ax)
        return ax

    def set_legend(ax):
        l = ax.legend(axespad=-0.05)
        plt.setp(l.get_texts(), fontsize=10)

    cmap_my = mpl.colors.ListedColormap(['g', 'b'])

    WS_false = []
    for i , j in enumerate(IR):
    
        if j >= 10.:    
            WS_it = 8.5
        else:       
            WS_it = 7.5
        
        WS_false.append(WS_it)
            
    ax = new_axes()
    ax.bar(Wdir[1:],WS_false, normed = True, opening = 0.8, edgecolor='white',nsector = 12,\
    cmap = cmap_my, bins = 2)

    set_legend(ax)
    
#-------------------------------------------------------------------------------

def SectorFrequencyPlot(T,RH,IR,Wdir):
    
# Function: SectorFrequencyPlot
#
# Purpose:  Plots sectorwise (12 sectors) icing intensity, relative humidity 
#           at T<0C with 0 icing intensity (Sublimation conditions) and Temperature
#           over 0 in a stacked plot. 
#           Inputs are Temperature, Relative humidity, Ice rate and wind
#           direction lists
#
# Note:     To be used with len(IR) = len(Wdir) - 1 ; The code uses Wdir[1:]
#
# Created:  31/01/2013    
        
    sect0 = []
    indT2 = []
    indT0 = []
    indRH70 = []
    indRH70_90 = []
    indRH90_95 = []
    indRH95 = []
    inddM10 = []
    inddM50 = []
    
    T = np.array(T)
    RH = np.array(RH)
    IR = np.array(IR)
    Wdir = np.array(Wdir)

    Dir_bins = np.arange(15,360,30)

    for p , k in enumerate(Wdir[1:]):
        
        if k <= 15. or  k > 345. :
            sect0.append(p)
                    
            if T[p] >= 2.:
                indT2.append(p)
                
            elif T[p] >= 0. and T[p] < 2. :
                indT0.append(p)
                
            elif RH[p] < 70. and T[p] < 0  and IR[p] < 10. :
                indRH70.append(p)
            
            elif RH[p] >= 70. and RH[p] < 90 and T[p] < 0  and IR[p] < 10. :
                indRH70_90.append(p)
                
            elif RH[p] >= 90. and RH[p] < 95 and T[p] < 0  and IR[p] < 10. :
                indRH90_95.append(p)
                
            elif RH[p] >= 95. and T[p] < 0  and IR[p] < 10. :  
                indRH95.append(p)
                
            elif IR[p] < 50:
                inddM10.append(p)
            
            else:
                inddM50.append(p)
    
    freq_vector = [len(indT2), len(indT0), len(indRH70), len(indRH70_90), \
                       len(indRH90_95), len(indRH95), len(inddM10), len(inddM50)]    
    bin_freq0 = float(len(sect0))
        
    probs0 = np.array(freq_vector)/bin_freq0*100.        
    pr = [probs0]

    for o , s in enumerate(Dir_bins[0:-1]):
        
        sect_v = []
        indT2_ = []
        indT0_ = []
        indRH70_ = []
        indRH70_90_ = []
        indRH90_95_ = []
        indRH95_ = []
        inddM10_ = []
        inddM50_ = []
                       
        for p , k in enumerate(Wdir[1:]):
        
            if k > Dir_bins[o] and  k <= Dir_bins[o + 1]:
                sect_v.append(o)
                
                if T[p] >= 2.:
                    indT2_.append(p)
                
                elif T[p] >= 0. and T[p] < 2. :
                    indT0_.append(p)
                    
                elif RH[p] < 70. and T[p] < 0  and IR[p] < 10. :
                    indRH70_.append(p)
                
                elif RH[p] >= 70. and RH[p] < 90 and T[p] < 0  and IR[p] < 10. :
                    indRH70_90_.append(p)
                    
                elif RH[p] >= 90. and RH[p] < 95 and T[p] < 0  and IR[p] < 10. :
                    indRH90_95_.append(p)
                    
                elif RH[p] >= 95. and T[p] < 0  and IR[p] < 10. :  
                    indRH95_.append(p)
                    
                elif IR[p] < 50 and IR[p] >= 10. :
                    inddM10_.append(p)
                
                else:
                    inddM50_.append(p)
             
        freq = [len(indT2_), len(indT0_), \
              len(indRH70_), len(indRH70_90_), \
              len(indRH90_95_), len(indRH95_), \
              len(inddM10_), len(inddM50_)]
              
        probs = np.array(freq)/float(len(sect_v))*100.                
        pr.append(probs) 
     
    ind = np.arange(len(Dir_bins))
    len_fr = np.arange(len(freq))

    bt = np.array([pr[u][7] for u in range(0,len(ind))])

    # Figure Plot

    color_vec = ['r','g','w','k','y','c','m']

    legend_label = ['$dM > 50 g/h$','$dM > 10 g/h$','$RH \geq 95\% $', \
                '$ 90\% \leq RH < 95\% $', \
                '$70\% \leq RH < 90\% $','$RH < 70\% $', '$0 \leq T < 2^{\circ}C$', \
                '$T \geq 2^{\circ}C$']
    plt.figure()
    Po = plt.subplot(111)
    Po.bar(ind, bt , width = 0.8 , color = 'b', label = legend_label[0])
                
    for i , j in enumerate(len_fr[0:-1]):
           
        Po.bar(ind, [pr[u][len(pr[1]) - 2 - j] for u in range(0,len(ind))], width = 0.8 , \
        bottom = bt , color = color_vec[i], label = legend_label[i + 1])
        
        bt = bt + np.array([pr[u][len(pr[1]) - 2 - j] for u in range(0,len(ind))])
               
    Po.set_ylabel('Frequency [%]')
    Po.set_xlabel('Sector')
    Po.set_xticks(ind + 0.8/2.)
    Po.set_xticklabels(['0', '30', '60', '90', '120', '150','180','210','240',\
              '270','300','330'])
    
    # Shrink current axis by 30%
    box = Po.get_position()
    Po.set_position([box.x0, box.y0, box.width * 0.7, box.height])
    
    # Put a legend to the right of the current axis
    Po.legend(loc='center left', bbox_to_anchor=(1, 0.5))          
    plt.show()
    
#-------------------------------------------------------------------------------    

def MonthFrequencyPlot(T,RH,IR,dates): 
    
# Function: MonthFrequencyPlot
#
# Purpose:  Plots monthly icing intensities, relative humidity 
#           at T<0C with 0 icing intensity (Sublimation conditions) and Temperature
#           over 0 in a stacked plot. 
#           Inputs are Temperature (list), Relative humidity(list), 
#           Ice rate(list) and dates (string in %Y-%m-%d %H:%M format) 
#
# Note:     To be used with len(IR) = len(Wdir) - 1 ; The code uses Wdir[1:]
#
# Created:  31/01/2013       

    m_id = []
    pr = []

    for i, j in enumerate(dates):
        m_id.append(dti.datetime.strptime(j, "%Y-%m-%d %H:%M").month)
        
    months = np.arange(1,13)    
        
    for o , s in enumerate(months):
        
        sect_v = []
        indT2_ = []
        indT0_ = []
        indRH70_ = []
        indRH70_90_ = []
        indRH90_95_ = []
        indRH95_ = []
        inddM10_ = []
        inddM50_ = []
                       
        for p , k in enumerate(m_id[1:]):
        
            if  k == s:
                sect_v.append(o)
                
                if T[p] >= 2.:
                    indT2_.append(p)
                
                elif T[p] >= 0. and T[p] < 2. :
                    indT0_.append(p)
                    
                elif RH[p] < 70. and T[p] < 0  and IR[p] < 10. :
                    indRH70_.append(p)
                
                elif RH[p] >= 70. and RH[p] < 90 and T[p] < 0  and IR[p] < 10. :
                    indRH70_90_.append(p)
                    
                elif RH[p] >= 90. and RH[p] < 95 and T[p] < 0  and IR[p] < 10. :
                    indRH90_95_.append(p)
                    
                elif RH[p] >= 95. and T[p] < 0  and IR[p] < 10. :  
                    indRH95_.append(p)
                    
                elif IR[p] < 50 and IR[p] >= 10. :
                    inddM10_.append(p)
                
                else:
                    inddM50_.append(p)
             
        freq = [len(indT2_), len(indT0_), \
              len(indRH70_), len(indRH70_90_), \
              len(indRH90_95_), len(indRH95_), \
              len(inddM10_), len(inddM50_)]
        
        if len(sect_v) == 0:
            probs = np.zeros(8)
        else:            
            probs = np.array(freq)/float(len(sect_v))*100.                
        
        pr.append(probs) 
     
    ind = np.arange(len(months))
    len_fr = np.arange(len(freq))
    
    bt = np.array([pr[u][7] for u in range(0,len(ind))])
    
    ## Figure Plot
    
    color_vec = ['r','g','w','k','y','c','m']
    
    legend_label = ['$dM > 50 g/h$','$dM > 10 g/h$','$RH \geq 95\% $', \
                    '$ 90\% \leq RH < 95\% $', \
                    '$70\% \leq RH < 90\% $','$RH < 70\% $', '$0 \leq T < 2^{\circ}C$', \
                    '$T \geq 2^{\circ}C$']
    plt.figure()    
    Po = plt.subplot(111)
    Po.bar(ind, bt , width = 0.8 , color = 'b', label = legend_label[0])
                    
    for i , j in enumerate(len_fr[0:-1]):
           
        Po.bar(ind, [pr[u][len(pr[1]) - 2 - j] for u in range(0,len(ind))], width = 0.8 , \
        bottom = bt , color = color_vec[i], label = legend_label[i + 1])
        
        bt = bt + np.array([pr[u][len(pr[1]) - 2 - j] for u in range(0,len(ind))])
               
    Po.set_ylabel('Frequency [%]')
    Po.set_xlabel('Month')
    Po.set_xticks(ind + 0.8/2.)
    Po.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun','Jul','Aug','Sep',\
              'Oct','Nov','Dec'])
    
    ## Shink current axis by 30%
    box = Po.get_position()
    Po.set_position([box.x0, box.y0, box.width * 0.7, box.height])
    #
    ## Put a legend to the right of the current axis
    Po.legend(loc='center left', bbox_to_anchor=(1, 0.5))          
    plt.show()
     
